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ABSTRACT 

Characterizing the surfaces of rocky exoplanets via their scattered light will be 
an essential challenge in investigating their habitability and the possible existence 
of life on their surfaces. We present a reconstruction method for fractional areas of 
different surface types from the colors of an Earth-like exoplanet. We create mock 
light curves for Earth without clouds using empirical data. These light curves 
are fitted to an isotropic scattering model consisting of four surface types: ocean, 
soil, snow, and vegetation. In an idealized situation where the photometric errors 
are only photon shot noise, we are able to reproduce the fractional areas of those 
components fairly well. The results offer some hope for detection of vegetation 
via the distinct spectral feature of photosynthesis on the Earth, known as the red 
edge. In our reconstruction method, Rayleigh scattering due to the atmosphere 
plays an important role, and for terrestrial exoplanets with an atmosphere similar 
to our Earth, it is possible to estimate the presence of oceans and an atmosphere 
simultaneously. 

Subject headings: astrobiology — Earth — scattering — techniques: photometric 
department of Physics, The University of Tokyo, Tokyo 113-0033, Japan 

2 Research Center for the Early Universe, Graduate School of Sciences, The University of Tokyo, Tokyo 
113-0033, Japan 

3 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA 

institute for the Physics and Mathematics of the Universe, The University of Tokyo, Kashiwa 277-8568, 
Japan 

5 Center of Climate System Research, The University of Tokyo, Kashiwa 277-8568, Japan 



-2- 



introduction 



A major milestone in exoplanet research will be the discovery of Earth-like planets. In- 
deed planets a few times as massive as the Earth, super- earths have already been discovered 
using the radial velocity method; for instance, M p sku of GJ581e is e stimated to ~ 2M m 



where M p is the planetary mass and i is the inclination of its orbit ( IMayor et al.ll2009l ). 
Microlensing and transit met hods are well-suited for the detection of such low-mass p lan- 
ets; MOA-2007-BLG-192-Lb faennett et al]|2008h and CoRoT-7b ffoueloz et all bopgh are 



reported to have masses of 3.3I^M ffi and 4.8 ± 0.8M®, respectively. The Kepler mission, 
launched on 2009 March 6, aims to detect a number of terrestrial planets with masses on the 
order of M e . Such low- mass planets are very likely to be rocky an d may have bodies of liquid 
water on their surfaces if they orbit in Habitable Zone (HZ, e.g., iKasting et al.lll993l 120031 ) 
of the primary star. Their discovery will definitely trigger ever more serious investigations 
of techniques to search for signatures of life, or biomarkers. 

The most conventional and extensively studied biomarkers are based on spectroscopic 
identification of molecular species in the planetary atmospheres, such as O2 and O3 (e.g., 
Leger et al.lll993l ; iDes Marais et al.ll2002l ). For the most favorable known exoplanetary sys- 
tems, detection of such atmospheric absorption features is already within reach during the 
transit or th e secondary eclipse; absorption features of Na, C, O, and H in a hot Jupiter 



HD209458b (Icharbonneau et al II2OO2I : Ividal-Madiar et al.lEooi 



2004 



12007 



Ballester et al. 



Swain et al. 



2008, 



200 



7) 



and H 2 0,CH 4 , CO, and C0 2 in HD189733b f lTinetti et all 
have been reported. Spectroscopy of non-transiting terrestrial planets orbiting in the HZ 
will likely require major and specialized observatories in space, but such facilities are being 
actively studied and evaluated by both NASA and ESA at present (e.g., 
Lawson et al.ll2009t iDarwin Mission Summary Status! 120071 ). 



Levine et al., 2009 



An even more ambitious and direct approach to the search for life may be possible via 
multi-band photometry of terrestrial exoplanets that is observationally less demanding. In 
particular, the light scat tered by th e surfa ce of a planet carries important information on 
properties of the surface. iFord et al.l ( 120011 ) computed for the first time the diurnal variation 
of scattered light by the Earth observed at a distance of 10 pc. They showed that the 
fractional variation of light curves of t he Earth is 10 % -2 0%, which indeed agrees with the 
result of Earthshine observations (e.g. 



Goode et al.l 120011 ) . 



More importantly, these model light curves exhibited an increase in brightness due to the 
presence of vegetation around at A =750 nm. This sharp increase of the albedo is known as 
the red edge. It is a fairly generic feature of vegetation on the Earth and is due to bio-pigments 
associat ed with photosynthesis. The red edge feature has been detected by Eart hshine obser- 
vations (j Woolf et al.1120021 Arnold et al.ll2002t iMontanes- Rodriguez et al.ll2006l ). It was sub- 
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sequen tly d iscussed as a possibl e bi omarker in extrasolar terrestr ial planets by ISeager et al. 



(2005) and iKiang et al.l (I2007al ]bf). iTinetti et al.l fj2006al Jbh and Montanes-Rodriguez et al. 



( 120061 ) performed more accurate simulations of the scattered radiation spectra of the Earth. 
They paid particular attention to the red edge, and discussed the detectability of its spectral 
feature. 



Palle et al.l ( 120081 ) focused on the determination of the rotation period of an extraso- 
lar terrestrial planet from the time variation of its scattered light. They concluded that 
the period can be reliably estimated in the presence of realistic partial cloud coverage of 
the planetary surface. They even found that the global circulation of the cloud pattern 
systematically modulates the estimate of the spin rotation period of the planet. 



More recently, ICowan et al.l (120091 ) approached the problem in a different and inter- 
esting manner. They performed a principal component analysis (PCA) on the real light 
curves of the Earth observed by the Extrasolar Planet Observation and Deep Impact Ex- 
tended Investigation (EPOXI) mission. They identified two major eigenspectra that roughly 
correspond to ocean and land on the Earth and extracted a rough distribution of these com- 
ponents as a function of longitude. This was the first attempt to solve the inverse problem 
and to extract the surface prop ertie s from the observed da ta in a model-independent fash- 



ion. I Williams &: Gaidosl ( 120081 ) and lOakley &: Cash! (120091 ) paid attention to the variation 
of scattered light according to the orbital motion of the planet. They pointed out that the 
reflectivity of ocean dramatically increases at a crescent phase a nd that the existence o f the 
liquid water may be observationally inferred through the effect. lOakley &: Cash! (120091 ) also 
suggested a possibility of reconstructing the land/ocean distribution along longitude using 
the gap between the reflectivity of ocean and that of land. 

In this paper, we develop still another method to estimate the fractional areas of dif- 
ferent surface types on extrasolar planets from multi-band photometry. We are particularly 
interested in terrestrial exoplanets in habitable zone and consider the detectability of vege- 
tation using the red-edge feature. Our attempt is to reproduce the scattered light curves as 
a sum of the four surface types, i.e., ocean, soil, vegetation, and snow. While PCA attempts 
to extract only orthogonal eigenspectra in a model-independent manner, the correspondence 
with real surface types is not clear. Indeed, we would like to answer an ambitious, but well- 
defined, scientific question: "If we discover an Earth-like exoplanet in the near future, to 
what extent can we reconstruct the surface information, in particular the presence of vegeta- 
tion observationally?". For that purpose, we intentionally and specifically adopt the major 
surface types of the Earth, including vegetation, into the analysis even though our method is 
thus inevitably model dependent. We discuss the feasibility of our method by applying it to 
mock scattered light curves of the cloudless Earth. We hope that our current results are rele- 
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vant for a future TPF mission (ILevine et al.ll2009f ) such as th e Occulting Ozone Observatory 
which aims at multi-band photometry of Earth-like planets (IKasdin et al.ll2010[ ). 



The rest of this paper is organized as follows. Section [2] describes the generation of mock 
light curves of the Earth. The methodology, assumptions, and results of our inversion model 
to reconstruct the fractional areas are presented in Section [3j We discuss the importance of 
the photometric band selection and comparison with PCA in SectionHl Finally we summarize 
the conclusion and discuss future directions of investigation in Section 



2. Mock scattered light curves of the Earth 

2.1. Configuration of the Star-planet-observer System 

This section describes our computation method for generating mock light curves that 
are analyzed in the following section. Figure [1] illustrates the geometric configuration of 
the star-planet system adopted in this paper. The vectors es, e^, and eo denote the unit 
vectors from the center of the planet (E) toward the host star (S), an arbitrary point on 
the surface (R), and the observer (O), respectively. A phase angle, a, is defined as an angle 
between es and eo, which is fixed to a = it /2 in this paper, except where noted otherwise. 
The solar zenith angle, 9q, at R is an angle between es and e^, and the zenith angle of the 
observation, 9±, is that between e^ and eo- The azimuthal angle between incident direction 
and the direction of observation is represented by cf) as shown in the right panel of Figure 
CD In our fiducial configuration, the planet is located at the equinox (except in Section I3.4f) 
and the observer is on the equatorial plane. 

The incident ray from the host star toward the planet first passes through the atmo- 
sphere, which may scatter the ray. A fraction of the light reaches the planetary surface and 
is scattered there, and then returns to space passing through the scattering atmosphere once 
again. The scattering properties of the atmosphere depend on its composition and overall 
optical thickness, while those of the surface vary according to the details such as compo- 
nents, coarseness, moisture, the shape of surface, and so on. The wavelength-dependence of 
the synthetic scattered light therefore carries information about the atmosphere as well as 
the surface features. The fact that the illuminated and visible part of the planet gradually 
changes in time due to the motion of the planet (orbiting around its host star and spinning 
about its own axis) allows observers to scan the planetary surface in a fashion similar to 
global remote sensing of the Earth observation satellites. Thus, the variation of planetary 
light curves in different photometric bands allows mapping the planetary surface to some 
extent. 
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spin axis 




Fig. 1. — Schematic illustration of our fiducial configuration of the planetary system. The 
star is located at the equinox and the observer is in the equatorial plane. The phase angle 
a denotes the planet-centric angle between the host star and the observer, which is fixed to 
a = 7r/2. The unit vectors eg and eo go from the center of the planet (E) toward the host 
star (S) and the observer (O), respectively. The vector e R is a unit vector normal to the 
surface plane at R. The solar zenith angle 9q, the zenith angle of the observation X , and 
the relative azimuthal angle are defined at each point R as the left panel indicates. As 
mentioned in Section |2T2| the illuminated pixels satisfy es ■ e^ > 0, and visible pixels satisfy 
e Q ■ e R > 0. 

2.2. Model and Assumptions 

To generate synthetic scattered light curves, we approximate the planetary sphere by 
polygons of 2°. 5 x 2°. 5 "pixels." The angles 6q, 0i, and </> are defined at the center of each 
pixel. In general, the radiance L(9 , 9\, <fi; A) (Wstr -1 m -2 /iin" 1 ) of each pixel is expressed 
as 

L(9 , 0i, 0; A) = F*(A) cos0 o /(0o, 9 U 0; A), (1) 

where F,(A) (Wm" 2 /im 1 ) is the incident flux at wavelength A, and f(9 , 9%, 0; A) (str -1 ) 
characterizes the scattering property of the surface and is conventionally referred to as the 
bidirectional reflectance distribution function (BRDF). The intensity of the total scattered 
light of a planet, /(A) (W str -1 /im -1 ), is obtained by integrating L over the illuminated and 
visible region of the planetary surface S: 

/(A) = / L(9 ,9 1 ,4>;X)cos9 1 dS 
Js 

= [ F»(A)/(0 o , 0i, & A) cos O cos O x dS 
Js 

= F*(\)R 2 p /"/(0 o ,0 1 ,0;A)cos0oCos0ids, (2) 

J s 
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where R p is the radius of the planet, ds = dS/R 2 is a normalized area element, and the 
illuminated and visible area S is defined as a set of pixels that satisfy e§ • > and 
e o " e R > simultaneously. Since the scattering property of a particular position on the 
surface gradually changes as the planet spins, 1(A) varies as a function of time. Equation 
02]) is a key formula in computing the synthetic scattered light curve. 

The evaluation of Equation (j2J) requires three inputs: R p , F*(A) and f(9 ,9i,(f>; A). We 
approximate the incident flux .F*(A) from the star by a black-body spectrum: 

FJX) = 2nkc2Rl 1 1 , (3) 

A 5 exp (he/ XksT^) — Id 2 '' 

where h is the Planck constant, c is the speed of light, is the Boltzmann constant, R* 
and T* are the radius and the surface temperature of the star, and d is the distance between 
the star and the planet. In this paper, we assume the Sun-Earth system, and therefore we 
set R p = R 9 = 6.4 x 10 6 m, T* = T = 5800 K, R* = R Q = 7.0 x 10 8 m, and d = 1 AU. 

We adopt a single scattering approximation for the interaction between the atmosphere 
and the underlying surface, i.e., the light is scattered once at most either in the atmosphere 
or at the surface and all multi-scattering processes are ignored. Then, the radiance of each 
pixel can be written as a linear combination of those associated with the surface and the 
atmosphere: 

L(6 , 0i, 0; A) = L atm (fl , 6 U 0; A) + L surf (0 O) Oi, <P; A), (4) 
or equivalently in terms of the BRDF (Equation ([T])), 

f(9 , 0i, 0; A) = / atm (0 o , 0i, 0; A) + f SUTi (6 , X , 0; A). (5) 
The specific models for / atm and / sur f are described below. 



2.2.1. / atm — Scattering by Atmosphere 

For scattering in the atmosphere represented by / a t m , we consider Rayleigh scattering 
alone, and ignore the effects of clouds, aerosols, and molecular absorption, which will be 
discussed elsewhere. We use a plane-parallel approximation locally for each pixel because 
our pixel size of 2.5° x 2.5° is sufficiently small that one can reasonably approximate the 
overall sphere by a polygon consisting of such pixels. While this approximation is not valid 
for pixels near the edge of the illuminated and visible area, the fractional area of the region 
is small and our results are not changed significantly. 
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Under the single-scattering approximation, the radiance from the atmosphere, L atm , is 
given by 



v atm 



or 



(9 Q ,9 l ,<j>;\) = F, / ^(e)exp 
Jo 

_ F«o;tt(e) 1 

cos Q\ 1/ cos 6*0 + 1/ cos Q\ 

w*(e) 



i i 

+ 



cos 6*o cos 0i 
1 



cxp 



COS 0i 

1 



/atm(0O,01,0; A) 



COS 0Q + COS 0i 



1 — exp < — r(A) 



COS 0o COS 01 

1 1 



(6) 



(7) 



COS 0o COS 01 

In the above expressions, u is a single scattering albedo that is the ratio of scattering efficiency 
to total light attenuation (both scattering and absorption). In what follows, we ignore the 
effect of absorption and assume u — 1. 



The phase function for Rayleigh scattering, \l/(0), is given by 

3 



*(6) 



167T 



1 + cos 2 6), 



where G is the angle between the incident and scattered directions: 

cos = cos 0o cos 0i + sin 0q sin 0i cos 0. 



(8) 



(9) 



The optical depth for Rayleigh scattering, r, is emp irically given for the atmosphere of the 
Earth (e.g., iFrohlich and Shaw lll98ol : hfoung Ill98f)l ): 



T 



0.00864 



P 



1013.25 mbar 



^-(3.916+0. 074A+0.05/A) ^ ^-4 



(10) 



where A is wavelength in units of 1 /im, and P is the pressure of the atmosphere at the base 
and we set P = 1013.25 mbar everywhere for simplicity. The fact that r is much less than 
unity in the visible and the near-infrared bands justifies the single-scattering approximation, 
that is adopted in this paper. 



2.2.2. / S urf — Scattering by Surface 

The radiance of the surface, L suri) is affected by atmospheric scattering as well, and is 
calculated by 

W0o,0i,0;A) = F*(A)cos0 o /surf(0O,01,0;A), (11) 

where 

/surf (0o, 0i, 0; A) = C atm (0o, 0i, 0; A)/ Sur f o(0o, 0i, 0; A), (12) 
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/surf, o is the BRDF of the surface itself, and C a t m describes attenuation due to atmospheric 
scattering before and after the ray reaches the surface: 

C atm = exp (-r(A) (-L- + -*—) ) . (13) 

L \ COS 6q COS d\ J J 

We consider the contributions of land, /land, cb an d ocean, / ocean , o, in Section [2 .2. 31 and Section 
I2.2.4| respectively. 



2.2.3. /i a nd, o — Scattering by Land 



Previous authors ( Ford et al. 2001 ;_ Tinetti et al. 2006al lbl; iMontanes-Rodrfguez et al~ 



2006 ; Palle et al. 2008 ; Oakley fc Cash 20091 ) first classified the pixels on the surface into dif- 



ferent surface types, and then computed the scattering due to each pixel according to its sur- 
face type. We do not apply any such classification, but rather use a parametrized BRDF for 
each pixel. More spe cifically, we adopt the mo del of MODerate resolution Imaging Spectro- 
radiometer (MODIS; ISalomonson et al.l (119891 )). an instrument on board the EOS/TERRA 
and AQUA satellites. The MODIS team processed the data using the Rossi-Li model for 
the BRDF, and assigned wavelength-dependent coefficients in this model to each observed 
pixel. The Rossi-Li model is one of sem i-empirical kernel-based BRDF models and has been 



applied widely (e.g. JLucht et al.l 120001 ) . It consists of three terms: 



/land(RL)(#0, <j>\ A) = /iso(A) + / v ol(A)K vo l(^0, #1, 0) + f SeO {X)K geo (B , &l, 0)- (14) 

The first term represents the isotropic component. The second term (the volume-scattering 
term) represents the effect of the finite thickness of the scattering body. The last term (the 
geometric-optical term) takes into account of the effect of shadow. The explicit expressions 
for K vo i and K geo are given in Appendix [A] 

In this paper, we adopt the three coefficients in the Rossi-Li model (/ iso , / vo i, / geo ) on each 
2.5° x 2.5° pixel from the dataset "snow-free gap-filled MODIS BRDF Model Parameters." 
This dataset is a spatially and temporally averaged product derived from the 0°.05 resolution 
BRDF/albedo data (v004 MCD43C1)Q. We select the BRDF parameters evaluated at the 
central position of each 2.5° x 2.5° pixel and employ the values to represent the pixel. 

The upper three panels of Figure [2] show the different behaviors of the surface brightness 
for each term in the Rossi-Li model. The lower-left panel plots /iand(RL)(^o> 4>'-i A) for 



1 The data-set is available through the MODIS web page 



http://modis.gsfc.nasa.gov/ 
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/iso = 0.236, / V oi = 0.114 and / geo = 0.027 (these values are obtained by averaging over 
the pixels corresponding to land). Note that /iand(RL) computed from Equation (}T4"]) is not 
positive definite, and becomes negative for very small #0 or 9\. Thus, we set /iand(RL) to 
wherever it becomes negative. The lower-right panel of Figure |2] shows the scattering surface 
brightness due to the ocean described below. 



2.2.4- /ocean, o — Scattering by Ocean 

The data set "snow-free gap-filled MODIS BRDF Model Parameters" does not assign 
parameters (/ iso , / vo i, / geo ) for ocean pixels nor for pixels around the polar region. Our 
current model regards those pixels as ocean. Thus, our model systematically underestimates 
the snow and/or ice areas around the poles. Those components actually vary significantly 
from season to season as well, which is not yet taken into account in our present model 
either. 

Since the scattering pattern by a liquid is very different from that by land, the above 
Rossi-Li model (Equation (fT4"|)) is not relevant for oceans. Instead, we adopt the ocean BRDF 



model of iNakajima fc Tanakal (119831 ) . We describe their model below. Further details are 



discussed in Appendix IBl 

Their model approximately computes the scattering of ocean by the sum of contributions 
from small facets (Appendix |B] and Figure [T8|) . The facets are characterized by the angle of 
their normal directions 9 n , and the scattering for each facet follows the simple Snell-Fresnel 
law. The distribution function p(9 n ) for the slope of those facets depends on the wind speed 
above the ocean (Equations (IB9P and (IBlOj) ) . Thus, their BRDF model is expressed as 

/ocean(#oA,0;A) = — — -p(9* n ) G(9 ,9 1 , 0, u 10 )r (9 ,9 1 , 91 0, m) , (15) 

4 cos 9 Q cos &i cos 6 1 * 

where #* indicates the direction of the wave slope responsible for the specular reflection, u w 
is the wind speed at 10 m height above the surface, and in = m(A) is the ratio of the refrac- 
tive index of ocean to atmosphere. The term r(9o,9i,(j),rh) stands for Fresnel's scattering 
coefficient, p(9 n ) is the density distribution function of the wave slope, and G(9 , 9\, 4>, Uio) 
represents the bidirectional shadowing effect. 

Our simulation adopts u w = 4 ms _1 , which corresponds to the spatial average of the 
wind speed at 10 m above the ocean, which is the average of monthly wind speeclE The 



2 the NOAA-ESRL Physical Sciences Division, Boulder Colorado. 



http: / / www.cdc.noaa.gov/ 



- 10 - 



color of the ocean is affected by the scattering of light inside the ocean, but we neglect this 
process for simplicity. We plan to incorporate it in future work. 

The surface brightness of ocean in this model is shown in the lower middle panel of 
Figure El The color contours clearly exhibit the remarkable increase in reflectance of the 
ocean near 8q = Q\ due to specular reflection. The strong specular reflection feature will be 
further discussed in Section [3~5l 
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Fig. 2. — Surface brightness of the sunlit hemisphere of the planet seen at a = tt/2 nor- 
malized by solar flux. The incident ray is coming from the right. The "isotropic" panel 
shows the surface brightness due to the isotropic term (assuming / iso = 1), the "volume" 
panel is due to the f vo \K vo \ term (assuming / vo i = 1), the "geometric" panel is due to the 
f geo K geo term (assuming / geo = 1), and the "Rossi-Li" panel is due to the combination of 
the three terms, /rl, assuming (/j SO , / vo i, / g eo) = (0.236, 0.114, 0.027), a set of averaged values 
over the entire land surfac e. The "ocean" panel is based on the ocean BRDF model by 
Nakajima &: Tanakal ( 119831 ) with a wind velocity 4 ms^ 1 . The "Rayleigh (SS)" panel is the 
surface brightness due to Rayleigh scattering of the atmosphere within a single-scattering 
(SS) and a flat atmospheric layer approximations, which are quantitatively not accurate at 
the edges, where either 9 and 6\ are close to zero. 



Adopting the above scattering model, we are now able to compute scattered light curves 
of a "Second Earth" (but without clouds). The planet rotates around its spin axis with a 
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period of 24 hr. The parameters of our mock observations are listed in Table [TJ 

We ignore the effect of the spin rotation during the exposure time £ exp . The orbital 
motion of the planet is completely ignored during the observation period n. We employ the 
set of MODIS photometric bands on which our data for the land BRDF are based. These 
bands are 0.459-0.479 ^m (band 1), 0.545-0.565 ^m (band 2), 0.620-0.670 /an (band 3), 
0.841-0.876 fjan (band 4), 1.230-1.250 /im (band 5), 1.628-1.652 /an (band 6), and 2.105- 
2.155 [im (band 7). In reality, we do not integrate over the band but simply calculate the 
scattering at the central wavelength of each band and multiply it by the band width. 

In practice, the spin rotation period of th e planet is unknow n a priori, and we must first 



determine it from the observed light curves. iPalle et al.l (120081 ) discussed how to infer the 



spin rotation period from the photometric variation and found that auto-correlation analysis 
of light curves can effectively determine the period. In this paper, therefore, we assume 
that the spin rotation period of the planet is precisely known, and fold the light curves 
accordingly. 



2.3. Results 



Figure [3] shows the simulated light curves in the seven bands. We consider a very 
idealized observational situation in which the light from the host star is completely blocked, 
and the photometric noise is due to the Poisson fluctuations in the observed photon counts 
from the planet alone: 



N ~ 840 ( — T — T ] ( — — 

VlO 15 Wstr-Vni- 1 / ^lOpc 




n 



14 days 




Table 1. Canonical Values of Parameters Assumed in Our Simulations of Photometric 

Light Curves. 



Parameter 


Symbol 


Value 


Planet-observer distance 


/ 


10 pc 


Effective diameter of telescope 


D 


2 m 


Exposure time 


texp 


1 hr 


Observation period 


n 


14 days 
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Fig. 3. — Light curves of the cloudless Earth at 0.459-0.479 /im (band 1), 0.545-0.565 /im 
(band 2), 0.620-0.670 /im (band 3), 0.841-0.876 /im (band 4), 1.230-1.250 /im (band 5), 
1.628-1.652 /im (band 6) and 2.105-2.155 /im (band 7). The error bars come from photon 
shot noise only. The bottom right panel shows the snapshots of the Earth at corresponding 
epochs. The ocean is painted in black and the land is in white. The observer is located on 
celestial equator and half of the projected planetary surface is illuminated. 



Here, the errors scale as ~ v N. In reality, however, there are certainly many sources of errors 
such as contamination by the host star, zodiacal light and detector noise; furthermore the flux 
from the planet will be attenuated through the telescope optical instrument and detectors. 
However for this preliminary investigation, we consider this idealized observational situation 
and leave a more realistic treatment of errors to future studies. 



Light curves at short wavelengths (especially in bands 1 and 2) in Figure [3] do not exhibit 
significant time variation because Rayleigh scattering by the atmosphere is dominant at 
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shorter wavelengths and dilutes the variations of surface features. The atmosphere becomes 
significantly more transparent at longer wavelengths and the variations of light curves due to 
the inhomogeneous surface become appreciable. A comparison of the light curves with the 
snapshots in Figure E] indicates that the three bumps at t ~ 1, 8, and 14 hr indeed correspond 
to the Eurasian, African, and American continents. The highest peak shows up when the 
Sahara desert emerges. A dip at t ~ 20 hr occurs when the illuminated and visible part is 
covered with ocean, since the reflectivity of land is higher than that of ocean, especially at 
longer wavelengths. Another dip at t ~ 13 hr occurs when the South American continent 
replaces the specular reflection poi nt. All these featu res ar e consistent with the l ight curve 
simulation of a cloudless Earth by iFord et al.l (l200lh and lOakley &: Cash! (120091 ). and are 
essential in extracting information concerning the surface features of a planet from its light 
curves. 
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Fig. 4. — Scattering properties of different types of surfaces. Cumulative fraction of scat- 
tered light is plotted against the corresponding fractional area defined as responsible for the 
scattered light. The lines at t =8, 20 hr are plotted with three references — ocean with 
wind velocity Wio = 4 ms _1 , Rayleigh scattering and Lambertian. Different colors represent 
different bands as in Figure [3] (band 1: blue, band 2: green, band 3: red, band 4: brown). 
The light scattered by the ocean is very localized. 

Figure H] illustrates the degree of localization of the source of the observed flux from the 
planetary surface. According to the geometry of the system that we adopt here, a quarter 
of the planetary surface area is illuminated and visible to the observer. We sort all the 
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pixels located in the illuminated and visible part of the surface according to the amount of 
scattered light per area, and compute the cumulative factional scattered light as a function 
of the corresponding fractional area. In Figure IH if all the pixels contribute to the light 
equally, the resulting plot would be a straight line reaching 100% at the fractional area of 
25%. In reality, however, the plot is slightly curved even in the case of isotropic scattering 
(Lambertian) due to the curvature of the global surface (Figure [2]). Since light from the 
ocean region comes from a very localized specular spot, the dotted line in Figure H] is very 
steep at small fractional area — the ocean spot equivalent to ~ 2% of the surface area is 
responsible for nearly the entire light from the ocean. 
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Fig. 5. — Variation of colors of Earth against time. Considered bands are 0.459-0.479 jira 
(band 1), 0.545-0.565 fim (band 2), 0.620-0.670 fim (band 3), and 0.841-0.876 /mi (band 
4). Here, "color" of band a and b is defined as C a b = — 2.5 log -j^ where I a and are the 
energy fluxes per wavelength in band a and band b. 

Essentially, the information that we utilize in reconstructing the fractional areas is color 
variations at different phases of the surface. To put it more clearly, we define the color 
between band a and band b as 

C ab = -2.51og^, (17) 

where I a and lb are the energy fluxes per wavelength in band a and band b (Equation 
(T5])). Figure [5] shows time variations of C12, C13, and C34. The trajectory on the C13 — C34 
plane is plotted in Figure [6] together with the typical colors of ocean, snow, vegetation, 
soil, and Rayleigh scattering which will be described in the next section and in Figure [71 
The trajectory indicates that the illuminated and visible part of the surface is almost fully 
covered by ocean at t = 20 hr. The other tips at t = 8 and 13 hr correspond to the African 
continent with the Sahara desert and the South American continent with the Amazon forest, 
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respectively. These color variations play a key role in estimating fractional areas of different 
surface types. 




Fig. 6. — Noiseless light curve trajectory on a color-color diagram assuming a black-body 
spectrum as the incident flux. The "color" of band a and b is defined as C & = —2.5 log 7 s - 

-* b 

where I a and lb are the energy fluxes per wavelength in band a and band b. The numbers 
denote the time (hours) corresponding to those in Figure EJ The symbols labeled by "ocean," 
"soil," "vegetation," "snow," and "Rayleigh (SS)" indicate the locations of the typical scat- 
tering colors of that surface type with atmospheric layer on it. Wavelength-dependent albedo 
of these five types are exhibited in Figure [7J 



3. Reconstruction of the fractional areas of surface types 

In this section, we describe our methodology for reconstructing the fractional areas of 
different surface types from multi-band photometry and present the results from the analysis 
of the mock light curves described in Section [2j 



3.1. Inversion Method 



Our basic strategy for reconstructing the surface features of a planet is to fit the pho- 
tometric light curves with an a priori model of planetary scattering. In doing so, we adopt 



two major simplifying assumptions. One is the Lambertian model (i.e., isotropic scattering) 
for the surface, and the other is that the surface consists solely of four different types (ocean, 
soil, vegetation, and snow) plus an atmosphere. 

The Lambertian surface is one of the simplest models of scattering with constant ra- 
diance against any geometry of incident and scattered rays. The BRDF for a Lambertian 
surface is simply given as 



where the subscript k denotes an index of the surface types and a^(A) is a wavelength- 
dependent albedo of the k-th surface type. 

The validity and limitation of our reconstruction method crucially depend on the number 
of different surface types that we consider. In practice, however, the limited information 
of color variations strongly restricts the number that can be uniquely determined by the 
analysis. Thus we consider only four types (k = ocean, soil, vegetation, and snow) that 
constitute the major components of the surface of the Earth. 

This approach inevitably limits the generality of our model. Nevertheless we think it 
reasonable for the present purposes for several reasons. First, our Earth is currently the only 
known example of a habitable planet, and it is not unreasonable to assume that at least some 
habitable exoplanets have similar surface properties. Second, these four types represent very 
different albedos (Figure [7]), which makes it easier to distinguish them in scattered light. 
Finally, the presence of ocean(s) on a planet is deeply related to the fundamental question 
of its habitability and the vegetation red edge can be regarded, if it exists at all, as a direct 
indication of the presence of life. 

Under the conditions described above, Equation ([2]) at a given epoch U reduces to 



where the integration is performed over Sk = Sk{ti) that is the area of the A;-th surface type 
in the illuminated and visible area at ti. Denoting each band by an index j (= 1, 2..., 7), we 
discretize Equation ( TT9l) as 





7T 



/(A) 




(19) 




(20) 
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Dj k — fiso jk > (21) 

TT 



A k {U) = ^ , (22) 



cos 6q cos 6ids 
cos 6> cos 9ids 

where A k is the geometrically-weighted fractional area of the k-th surface type that we want 
to estimate, and Dj k is normally referred to as the "design matrix" by statisticians. 

In Equation (|20|) . the incident flux F*j is calculated from the intensity of the host star 
and the distance between the star and the planet d, which is obtained once the orbit of the 
planet is determined. The integral J s cos0ocos0ids depends on the phase angle a alone (in 
our current configuration of the phase angle a = tt/2, J g cos0 cosOids = 2/3). 

In order to estimate the weighted fractional area A k (ti), we need to determine the design 
matrix, Dj k from the albedo cij k . We define the effective albedo of the k-th surface type, 
Q efffc(A), by setting it equal to the actual reflectivity when the whole surface is covered by 
that surface type. More specifically, it is computed as 



/model fc (00,01,0; A) COS 6» COsOids 

a cSk (\)=n^ ? , (23) 

cos^o cosOids 



where / mo dei k is the model BRDF of the k-th surface type and the integrations are performed 
over the illuminated and visible area which, again, depends on the phase angle a alone. Then 
we define cij k of the j'th band from its central wavelength, A-,-, as 

a jk = a cSk (Xj). (24) 

The model BRDF / mo dei k in Equation (1231) is calculated from 

/model k (00, 01, 0; A) = /atm(0O, 01, 0; A) + C atm (0O, 01, 05 A) / su rf ,0 k (00 , 01 , 05 A), (25) 

where / atm and C atm are given by Equations ([7]) and (fl3l) . respectively. In order to deter- 
mine /surf ,o k of the three land components (k = soil, vegetation and snow), we first assume 
th at they are Lamberti an with scattering spectra given by the ASTER spectral library 
|f| tealdridge et al. 2008 ). Specifically, we adopt "Brown to dark brown sand (Entisol)," 



"Grass," and "Fine Snow" for soil, vegetation, and snow, respectively. Since the ASTER 



3 http://speclib.jpl.nasa.gov/ 
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spectral library offers data at discrete wavelengths, we linearly interpolate the data to ob- 
tain the /surf, OA: suitable for our MODIS bands. For / sur f 5 oA: of ocean, we use Equation ffT5|) . 
Figure [7] displays the effective albedo at P = 1013.25 mbar for ocean, soil, vegetation, and 
snow as a function of wavelength A. The dashed lines are the albedo curves in the absence 
of an atmosphere. The black solid line is the effective albedo due to atmospheric Rayleigh 
scattering. The primary effect of Rayleigh scattering is to add a very blue continuum to 
every pixel's contribution to the total light, but in our single-scattering approximation it 
also reduces the amount of light reaching the surface. In Figure [TJ this is most clearly visible 
in the blue spectral region for the snow component. 



(1)(2)(3) (4) (5) (6) (7) 




wavelength [ixm] 

Fig. 7. — Wavelength-dependent effective albedos of ocean (blue), soil (magenta), vegetation 
(green), snow (cyan), and atmosphere with Rayleigh scattering alone (black). The solid 
lines show the effective albedo with Earth-like atmosphere, while the dashed lines show the 
effective albedo without an atmosphere. Shaded regions correspond to the MODIS bands. 
The numbers at the top are the labels of the different photometric bands. 



3.2. Fitting 

Now we apply the above methodology to the mock light curves (Section [2]). We deter- 
mine the best-fit values of A, (tj) in Equation ([201 by chi-square fitting: 

2/, n _ {tobsjjti) - F*jRp J s cos fl cos Oids J2k D jk('r)A k (t l )} 2 

The fit is independently performed for the light curves at each epoch t{. In order to quantify 
the errors of fitted values, we generate 100 realizations by adding a Poissonian error with 
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rms of a, to the simulated data, and compute the average and the variance of the best-fit 
values. We consider the photon shot noise alone in crj, ignoring the other statistical errors 
of the phase angle a, star-planet distance d, and planet-observer distance I among others. 
Thus our model is admittedly very idealized but indicates what one can learn from possible 
future data in principle. 

Equation (|26|) is essentially a linear inverse problem. So as to ensure the positivity 
in each element of Aj, however, we replace Aj by B? and search for the best-fit B~ with 



a nonlinear fitting method (the Levenberg-Marquardt method, e.g., iPress et al.lll992l ). We 



also made sure that another independent fitting method, which i s based on a conditional 



least square method and does not use the above replacement (e.g.. lMenkdll989l ). gives very 
similar results. 

Our model does not assume the condition ^2 k=1 A k = 1. We could impose that con- 
straint in principle, but it would somewhat restrict the applicability of the model. After all, 
the actual planetary surface does not exactly consist of only four Lambert surface types as 
we assumed here, and components other than the four types might contribute. In addition, 
the sum of the estimated area can easily deviate from unity due to a variety of other ef- 
fects including the anisotropy of the scattering and the diversity of the detail features of the 
reflection spectra (e.g., the sharpness of the red edge, the slope of the soil spectrum etc.), 
Therefore, we fit the mock data without any constraints other than A^iti) > 0. 



3.3. Estimation of the Weighted Fractional Areas 

In reality, the atmospheric optical depth due to Rayleigh scattering r will not be known 
in advance and should be fitted simultaneously from the data. In this subsection, however, 
we simply set the optical depth to the standard value r : 

r (A) = 0.00864A" 4 , (27) 

which is an approximation form of Equation f lTUj) with P =1013.25 mbar. The effect of r 
will be discussed in Section 13.51 

The top panel of Figure [8] shows the % 2 /d.o.f. (d.o.f. = 5 — 4 = 1) of the fitting at each 
epoch. The second panel from the top presents the result of estimating weighted fractional 
areas A k {ti) of ocean and land (= soil + vegetation + snow) using the light curves in bands 
1, 2, 3, 4, and 5. The third panel from the top panel displays the fractional areas of the 
three land types separately. The bottom panel illustrates the corresponding snapshot of 
the Earth toward the observer. The symbols indicate the average of the best-fit values 
from the 100 realizations with quoted error bars representing the rms among them. For 
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reference, we plot in dashed curves the weighted fractional areas based on the One-Minute 
Land Ecosystem Classification Product, which is generated from the official MODIS land 
ecosystem classification dataset classifying the surface of the Earth into 16 classes. Among 
the 16 classes, we regard "water" as ocean, "snow and ice" as snow, "open shrubland", 
"permanent wetlands", "urban and built-up", and "barren or sparsely vegetated" as soil, 
and others as vegetation, as shown in Table [2j 

The comparison of our best-fit values and dotted lines indicates that the present method 
works fairly well. Given the relatively crude approximation of the isotropic scattering and the 
assumption of only four surface types incorporated in the analysis, it is perhaps surprisingly 
successful. The presence of ocean, soil and vegetation is recovered, and the variation of their 
weighted fractional areas follows ground truth at least qualitatively. The bumps of the ocean 
at t ~ 11 and 20 hr correspond to the Atlantic Ocean and the Pacific Ocean, respectively. 
The peaks in the soil and vegetation curves indeed correspond to the Sahara desert and the 
Amazon forest, respectively. Moreover, the weighted fractional area of snow is consistent 
with zero, in agreement with the fact that we have adopted the snow-free land BRDF data 
in the light curve simulation. 

The quantitative discrepancy between the expected and model fractional areas of soil 
and of vegetation probably comes from the diversity of land scattering properties that is 
actually not well represented by the four Lambertian types shown in Figure [7J Note that 
our assignment of the 16 surface classification into the four types (Table [2]) is not unique, 
and therefore the dashed lines should be regarded as a plausible but not unique reference. 



Given the fact that scattered light from oceans comes almost exclusively from a very 
small region of specular reflection, it is perhaps puzzling that our inversion method, based 
on an assumption of Lambertian scattering, can estimate the fractional area of ocean rea- 
sonably well. This can be understood as follows: except for the specular reflection spot, 
the scattered light from the ocean surface is negligible (Figured]). Rayleigh scattering in 
the atmosphere adds a uniform and very blue continuum to the observed flux from every 
part of the illuminated and visible part of the planetary surface, including areas of ocean 
away from the specular reflection point. For land areas, the contribution to the scattered 
light is dominated by direct scattering from the surface, especially in the redder bands, so 
the effect of Rayleigh scattering is merely to produce a slightly bluer overall color (Fig. [7J). 
For ocean areas outside the specular reflection spot, however, the contribution to the total 
scattered light comes primarily from the Rayleigh scattering. Furthermore, the resulting 
scattered light by the atmosphere is fairly diffuse and comes from the entire illuminated and 
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Table 2: The International Geosphere-Biosphere Programme(IGBP) classification which is 
generated from the official MODIS land ecosystem classification data-set(MOD12Ql). The 
fourth column is the merged classification we adopt for dashed lines in Figure [H 
No. IGBP Classification Area(%) Our Classification Ocean/Land 






water 


71.40 


ocean 


ocean 


1 


evergreen needleleaf forest 


1.13 


vej 


fetation 


land 


2 


evergreen broadleaf forest 


2.87 


ve; 


fetation 


land 


3 


deciduous needleleaf forest 


0.18 


vej 


fetation 


land 


4 


deciduous broadleaf forest 


0.46 


vej 


fetation 


land 


5 


mixed forest 


1.34 


vej 


fetation 


land 


6 


closed shrubland 


0.16 


vej 


fetation 


land 


7 


open shrubland 


5.22 




soil 


land 


8 


woody savannas 


2.15 


vej 


fetation 


land 


9 


savannas 


1.99 


vej 


fetation 


land 


10 


grasslands 


2.65 


vej 


fetation 


land 


11 


permanent wetlands 


0.06 




soil 


land 


12 


croplands 


2.56 


ve; 


fetation 


land 


13 


urban and built-up 


0.14 




soil 


land 


14 


cropland/natural vegetation mosaic 


0.60 


ve; 


fetation 


land 


15 


snow and ice 


3.16 




snow 


land 


16 


barren or sparsely vegetated 


3.92 




soil 


land 


17 


unclassified 


0.00 
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Fig. 8. — Reconstructed fractional areas A^tj) for four surface types from the simulated 
light curves in five bands (bands 1-5). The top panel shows the value of reduced \ 2 
(=X 2 /dof where degree of freedom (dof) is 5 — 4 = 1) for each epoch. The upper mid- 
dle panel displays the results of estimating weighted fractional areas of ocean (blue), land 
(=soil+vegetation+snow; brown), and the total of them (red). The lower middle panel 
displays those of soil (magenta), vegetation (green), and snow (cyan). The dashed lines in 
those two panels show the weighted fractional areas derived from the MODIS land ecosystem 
classification dataset. The quoted error bars indicate the variance of the best-fit values from 
100 realizations. The bottom panel depicts the snapshots of the Earth at the corresponding 
epochs where the ocean is painted in gray and the land in white. 
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visible area. Thus, in order to account for the observed amount of Rayleigh scattering color 
in the total, the fit must assign a sufficiently large area of dark (very low albedo) area to 
the surface. Ocean is the only low albedo component available to the fit among the four 
components and thus is automatically assigned to produce enough Rayleigh scattering color 
without producing too much land scattering color. In other words, in our fit (which makes 
no use of spatial or time-domain information), the Rayleigh scattering component without 
any corresponding land component is a proxy for the ocean. The dark ocean surface is de- 
tectable in the fit to the colors primarily due to the additional Rayleigh scattering continuum 
supplied by the atmosphere above it. 



3.4. Time-integrated Spectra 

Since our fiducial observational condition (Table [I]) is fairly idealized, the quoted errors 
are small even in the case of time-dependent analysis. In an actual missions, however, time- 
intergrated spectra will probably be the first realistic goal. The spectrum from the total 
integration time of the mock observation (2 weeks) is displayed in Figure [9j This result is 
obtained by adding up the photon counts for each band over two weeks and then normalizing 
them with respect to the scattered intensity of a lossless Lambert sphere. The error bars 
come from the shot noise of the total photon count in each band. Note that the red edge 
produces a sharp increase at A ~ 800 nm. 
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Fig. 9. — Spectrum of the cloudless Earth from the total integration time (2 weeks). The 
i/- axis is normalized by the scattered light by a lossless Lambert sphere at full phase. 



We present in Figure [10] the result of the fit to this time-averaged spectrum with our 
fiducial configuration. The fitting model is now the time-averaged version of Equation (1201) 



-24- 



CO 
CD 



CO 

c 
o 

o 

CO 



100 
80 
60 
40 
20 



best-fit h 
actual weighted fraction 
actual fraction c 



ocean snow vege soil land 



Fig. 10. — Time- averaged spectrum from the whole integration time (Figure [9]) is decomposed 
into ocean, snow, vegetation, and soil with Equation ( 120]) . The land is the summation of 
snow, vegetation, and soil. The green boxes are the weighted fractional areas based on our 
merged classification (Table |2]) and the blue boxes show the non-weighted actual fractions. 



and thus the estimated values correspond to A^(tj) averaged over the full integration time. 
The blue boxes represent the actual fractional areas, while the green boxes represent the 
fractional areas based on our merged classification (Table |2]) but correspondingly weighted 
to be compared with the estimated value. Under our fiducial configuration, the fractional 
areas are reasonably reproduced even without time resolution (and hence without spatial 
resolution) . 

We also repeat the same fitting procedure using mock light curves with different ge- 
ometries, and summarize the results in Table [3j The values in the "reference" lines indicate 
weighted and averaged fractional areas visible for each observer. The bottom line shows the 
actual (non-weighted) fraction of the four surface components. Of course, the correspon- 
dence of the actual non-weighted fractional areas and the weighted fractional areas is highly 
dependent on the geometry. In the case of the Earth, the weighted fractions are roughly 
equivalent to the actual fractions in our fiducial configuration i.e. if the Earth is at equinox 
and the observer is on the equatorial plane, due to the fact that areas of ocean, soil, and 
vegetation are not extremely inhomogeneous but distributed well along the latitude. On the 
other hand, land areas would be preferentially seen as viewed by an observer located above 
the northern hemisphere. Although such geometric limitation is inevitable, the estimated 
values recover the weighted fractional areas for each geometry reasonably well. This is very 
encouraging for future missions. 
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Table 3. Estimated values of weighted fractional areas from time-averaged spectra with 
different geometries. The "north 45°" and "south 45°" assumes the Earth at equinox seen 

with +45° and —45° inclination from the equatorial plane, respectively. The "summer 
(winter) solstice" assumes the Earth at summer (winter) solstice and the observer on the 
intersection of the equatorial plane and orbital plane. The values in "reference" lines are 

based on our merged classification shown in Table [2] which are weighted and averaged 
according to each geometry. The "actual (non- weighted)" is the no n- weighted fractional 
areas of each components based on our merged classification. 



Geometry 




ocean (%) 


land (%) 


soil (%) 


vegetation (%) 


snow (%) 


fiducial 


estimated 


75.5 


21.5 


12.4 


8.3 


0.7 




reference 


74.0 


26.0 


8.8 


16.8 


0.3 


north 45° 


estimated 


66.5 


29.4 


17.5 


10.1 


1.7 




reference 


56.9 


43.1 


17.2 


25.1 


0.8 


south 45° 


estimated 


86.6 


10.7 


4.3 


5.6 


0.7 




reference 


85.6 


14.3 


4.0 


6.9 


3.4 


summer solstice 


estimated 


73.1 


26.5 


16.3 


9.4 


0.8 




reference 


68.0 


32.3 


11.5 


20.3 


0.03 


winter solstice 


estimated 


83.7 


16.6 


9.3 


7.3 


0.1 




reference 


79.2 


20.8 


6.3 


13.1 


1.4 


actual (non-weij 


dited) 


71.4 


28.6 


16.1 


9.3 


3.2 
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3.5. Estimation of the Optical Depth of Atmosphere 

So far we have assumed that the optical depth r is indeed equivalent to the input value 
To. In reality, however, the value of r would not be known a priori, and thus r should be 
regarded as one of the fitting parameters. The use of an incorrect value of r would degrade 
the fit. We investigate this issue by repeating the fit with different input optical depths: 

Tl 

T- fit = — r (ra = 0,1, 2,. ..,20). (28) 

The matrix Djk depends on the value of r (Equations (JTj), (EES]) , ( 1231) . and (1251)). Assuming 
that the optical depth is constant over the planetary surface, we sum up Equation (1261) over 
different epochs according to 

24 5 \ ^obsi(^i) - F*jRlf s cos 9 cos Oids ^2D jk (r)Ak(ti) > 

*V) = EE- m k ~ L < <»> 

and search for the best-fit value of r that minimizes Equation fl29|) among the different values 
(Equation ([28])). 



1.4 rr 




4 8 12 16 20 24 
time [hour] 

Fig. 11. — Dependence of reconstructed fractional areas on the input value of r^ t . The 
painted area exhibits the range of the estimated values when t^/tq varies from 0.5 to 1.5. 
Simulated light curves in bands 1 to 5 are used. The dashed lines show the weighted fractional 
areas derived from the MODIS land ecosystem classification dataset, like those in Figure |HJ 

It would be instructive to see first how the reconstructed fractional areas are sensitive 
to the value of rgt . Figure [TT] shows a case for which the simulated light curves with r = tq 
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are fit using input values of ret varied from 0.5ro to I.Sto- Since different values of rat 
mainly modify the effective albedo of ocean through Rayleigh scattering, estimates of ocean 
fraction are sensitive to rgt! ocean area is overestimated (underestimated) for small (large) 
ret- However, the estimate of the fractional areas for the other three surface types is fairly 
robust, indicating that these fractional areas are determined mainly at longer wavelengths 
where the value of r does not make any significant difference. 

Figure [T2] displays the histogram of the best-fit values of r fit (Equation ([25]) ) which 
minimize Equation (1291) . The results show a broad peak around r^ t ~ 1.2r , but cases of 
ret < 0.9t or rat > 1-6t are relatively rare. Given the crude approximations adopted in our 
reconstruction method, these results are also encouraging. 




Fig. 12. — Histogram of best-fit values of r& t from 100 realizations of simulated light curves 
in bands 1 to 5. 

Thus we have found that our inversion method can recover the presence of ocean and 
atmosphere simultaneously for a cloudless Earth-like planet. It may be also instructive to 
consider a hypothetical ocean planet similar to Earth but without atmosphere at all. In 
this case the ocean contributes a tiny fraction of the total scattered light, and also a small 
fraction in terms of area because of the specular reflection. 

In order to see this more quantitatively, we create light curves for the Earth without an 
atmosphere (i.e., Rayleigh scattering is neglected), and repeat the same analysis. The best- 
fit result with r fit = in Equation (J2"fi]) is shown in Figure [T3J As expected, the estimate 
of the fractional areas of ocean is very unstable and unreliable. It is interesting to note, 
however, that one can still reconstruct the fractional areas of soil and vegetation fairly well. 
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Fig. 13. — Same as Figure but without an atmosphere. 

4. Discussion 

4.1. Band Selection 

Our inversion method relies entirely on the difference of the wavelength-dependence 
of the scattering spectra among ocean, soil, vegetation, and snow, or their colors in short. 
Therefore the selection of the observed bands is crucial. 

First we repeat the same analysis performed in Section [3] but with four bands. In the 
case of five bands, we can fit up to five unknowns; we selected fractional areas for the four 
types and r as the five fitting parameters. In the case of 4 bands, however, we cannot fit 
t simultaneously and thus fix r = r . The result is shown in Figure HH The left and right 
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panels correspond to the bluer bands (j = 1, 2, 3, and 4) and redder bands (j = 4, 5, 6, and 
7), respectively. The difference can be easily understood; the bluer bands are more sensitive 
to the effect of Rayleigh scattering and the fractional area of ocean is recovered fairly well. 
Also the red edge feature between bands 3 and 4 still carries the vegetation signature, while 
the fractional area of soil becomes more uncertain because it is brighter in the redder bands. 
The result with the redder bands shows consistently opposite generic features; the lack of 
the blue bands and the red edge makes it difficult to detect the signature of ocean and 
vegetation, respectively, while the soil is more easily reconstructed. 
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Fig. 14. — Reconstructed fractional areas using simulated light curves in two different sets 
of four bands; Left: bluer four bands (j = 1, 2, 3, and 4), Right: redder four bands (j = 4, 
5, 6, and 7). 

A closer look at Figure [HI indicates some degeneracy between reflection spectra of 
selected surface types (Figure [7]). This is more clearly illustrated in Figure [151 where we use 
three bands (j = 2, 3, and 4) only. In this case we cannot determine four surface types, 
and we choose ocean, soil, and vegetation in the left panel, and ocean, vegetation, and snow 
in the right panel, while we neglect the remaining surface type from the fit. The result 
naturally is degraded compared with the 4-band and 5-band cases. The vegetation signature 
is still there because we chose bands 3 and 4 that bracket the red edge, and the fractional 
areas of soil and snow compensate for each other. The fact that the spectra of snow and 
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soil are similar in shorter wavelength bands (except for their amplitude) partly explains the 
behavior of the lower panels in Figure [T5j While this similarity may come largely from our 
single scattering approximation as described in Section 13. 1\ we obtained a similar result even 
when we use the snow effective albedo of r = 0. Therefore the degeneracy between snow 
and soil is fairly generic as long as we use the short wavelength bands in the reconstruction. 
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Fig. 15. — Reconstructed fractional areas using simulated light curves in three bands; Left: 
ocean, soil, and vegetation are considered. Right: ocean, vegetation, and snow are consid- 
ered. 



4.2. Principal Component Analysis (PCA) 



Cowan et al.l ( 120091 ) performed a PCA of the multi-band light curves of the Earth 



as observed by the EPOXI mission. They extracted two major eigenspectra in a model- 
independent manner. W e also performed PC A on the same data-set, and obtained results 
consistent with those of ICowan et al.l (120091 ) . We then performed PCA of our mock light 
curves and again extracted two major eigenspectra; these are shown in Figure | 



In order to compare our current methodology with PCA, we decompose these two eigen- 
spectra into a combination of the effective albedo of the four surface types with r = r 
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(solid lines in Figure [7]). We find that the first eigenspectrum roughly corresponds to 
(soil+vegetation-ocean), and the second one roughly corresponds to (vegetation- (soil+snow+ocean)) 
as displayed in Figure [T71 

This result indicates that the extracted eigenspectra do not necessarily correspond to 
any single surface type. This is not surprising, of course. While PCA extracts orthogonal 
eigenspectra by definition, the wavelength dependence of albedos of real materials is not 
orthogonal in general. Moreover, they are likely to be degenerate in PCA because the 
fractional areas are complementary; for example, the fractional area of ocean decreases when 
fractional area of land increases. In other words, the time dependences of these components 
are necessarily correlated. Therefore, it seems natural that the first eigenspectrum is roughly 
(land-ocean). Our method is quite model dependent, but it allows decomposition of the light 
curves into physical components that can be interpreted in a simple way. While the model 
independence of PCA is a great advantage, the final interpretation is not straightforward. 
Further comparison with PCA is beyond the scope of this paper, but clearly these two 
methodologies are very complementary. 
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Fig. 16.— Left: the eigenspectra extracted by PCA of our fiducial mock light curves 
displayed in Figure El The contribution of the first eigenspectrum (solid) is 94.3 % and that 
of the second one (dashed) is 5.7 %. Right: the time variation of the first eigenspectrum 
(solid) and the second eigenspectrum (dashed). 



5. Summary 

In this paper we have presented a method to reconstruct the fractional areas of different 
surface components from photometric colors of Earth-like exoplanets without clouds. For 
this purpose we first created mock light curves in seven photometric bands from the observed 
data of the Earth, but neglected the clouds. The light curves are fitted to isotropic scat- 
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Fig. 17. — Left: decomposition of the first eigenspectrum of Figure ITBl where the coefficients 
are ocean, -2.19; snow, 0.03; vegetation, 0.86; soil, 1.07. Effective albedos with r = r are 
used. Right: decomposition of the second eigenspectrum of Figure [TBI where the coefficients 
are ocean, -1.83; snow, -0.58; vegetation, 2.36; soil, -1.54. Effective albedos with r = r are 
used. 



tering models consisting of four surface types: ocean, soil, snow, and vegetation. In a very 
idealized situation where the data are obtained using a noiseless 2 m telescope and multiple 
integrations of 1 hr each, we find that our method is able to reproduce the fractional areas of 
those components fairly well. In particular, Figures IHl fTTl IT3| [Ml and [TBI show quantitatively 
that the presence of vegetation can be recovered using the color information via its red edge 
feature. Although we fit the fractional areas of each component independently at each time 
in the light curves, the time variation is event ually translated into t he distribution along the 



longitude with the methodology described in lCowan fc Agoll (120081 ). 



We also find that Rayleigh scattering due to the atmosphere plays a key role in esti- 
mating reliably the fractional areas of ocean. Indeed, without an atmosphere, our method 
based on the isotropic scattering assumption cannot properly estimate the real fractional 
area of ocean because of the strongly anisotropic nature of its specular reflection. On the 
other hand, for terrestrial exoplanets with atmosphere similar to our Earth, we may be able 
to estimate the presence of ocean and atmosphere simultaneously if the effect of clouds is 
safely neglected. 

Our methodology described in this paper is based on admittedly several very idealized 
assumptions and simplifications. There are a variety of issues that we have to address and 
improve in future work, some of which are listed below. 

First, one of the most serious omissions in the present modeling is the absence of clouds. 
Clouds provide additional time variation in the light curves, which is not directly related to 
the property of the planetary surface. Clouds greatly increase the reflectivity in visible and 
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near-infrared bands. We may treat clouds as an additional component in our current model. 
While the time and spatial variabilities of clouds are not easy to model, combining data from 
multiple rotation periods may allow us to separate variations due to surface features from 
those caused by clouds. We will discuss these effects of clouds more quantitatively in future 
work. 

Second, the universality of the scattering spectra of surface types on Earth is not clear, 
maybe not even likely. Indeed this could be regarded as both the strongest and the weakest 
point in our methodology; we can recover surface information, including the presence of 
vegetation from the color variations alone because our spectral template is derived from 
that of the Earth. Nevertheless, as stressed in the introduction, it is an important first 
attempt to see if we could infer the presence of vegetation observationally for exoplanets 
similar to the Earth even in the most favorable and idealized circumstances. The answer to 
the question seems promising and the next step is to generalize the result for a wider range 
of possibilities. The scattering properties of ocean and pure snow may be universal because 
the properties of H 2 would not be different on different planets, but they can be somewhat 
different depending on impurities, for example. The generality of the red edge depends on 
the photosynthesis system of exo-vegetation and t he wavelength of that "edge" is likely to 



shift a ccording to the spectral type of the host star (jWolstencroft fc Ravenll2002l ; iKiang et al. 



2007cM . Also, the scattering properties of soil may vary depending on the details of the soil's 



composition, and even the gravity may affect the granularity which changes the scattering 
properties. 

Finally we can explore other geometrical configurations of a star-planet system includ- 
ing the effects of orbital inclination, planet obliquity, and seasonal variations. It will be 
interesting to apply our methodology to time-series remote sensing data-sets for the Earth 
as well as to models of the Earth in the distant past which take into account continental 
drift and the evolutionary history of vegetation. 
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A. BRDF for land 

We create mock light curves adopting the Rossi-Li model (Equation ( JT4l) ) for the land 
BRDF. The specific expressions for the volume-scattering term, K vo \, and the geometric- 
op tical term, are given here for completeness. The derivation of the two terms is found 



in 



Wanner et al.l ( 119951 ). The volume-scattering term is 



v fa a a\ (tt/2-Q + sin£ vr 

K vol (6 ,8 U (j)) = — -, (Al) 

COS Uo + COS 9\ 4 

where 

cos £ = cos 6*o cos 9\ + sin O sin 9i cos <fi. (A2) 
The geometric-optical term is: 

> - sec9[ + -(1 + cos O sec 0q sec 0i, ( A 3) 



where 



O(9 ,9 1 ,4>) = -(t-sintcost)(sec0 o + sec0i), (A4) 
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r 



b 

-tan0i). (A9) 



r 



This kernel assumes a sparse ensemble of surface objects which throws shadows on the 
Lambertian background. The objects are approximated as spheroids with width 2r and 
vertical length 2b. The height of the center of the spheroids is denoted by h. The function 
O(0o, 0i, 0) is the overlap area between the solar shadows and the shade of view. For MODIS 
processing h/b = 2 and b/r = 1 are assumed (i.e., the spherical crowns are separated from 
the ground by their radii). Thus, we adopted these values. 
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B. BRDF for ocean 

The BRDF model for ocean is very com plicated, and we summarize the key expressions 
of the model by iNakajima fc Tanakal (119831 ) that we adopt in the present work. 



The scattering of solar radiation by a flat ocean follows the simple Snell-Fresnel law. 
Using the pair {9, <$} to represent polar coordinates relative to the normal of the averaged 
surface plane (Figure [T8]) . the relation between the incident radiance L m {6' ', 0') and the scat- 
tered radiance I OVLt (9,d>) is written as 



'out i 



R(9,9', 



d cos 9'd(f)'R(9,9',0 
r(9, m)5(cos 9 — cos 9')5(<f)' — <p — tt) 



(Bl) 
(B2) 



where r(#,m) is the Fresnel scattering coefficient and fh the is the ratio of the reflective 
indices of air and ocean. When the wind above ocean is taken into account, however, the 
slope of scattering surface varies both temporally and spatially. 



emergence 




Fig. 18. — Configuration of incident and scattered rays with respect to a wave facet whose 
normal vector n is {9 n , <ft n } in the polar coordinates with respect to the average surface 
plane. 

Thus, the normal direction of the facet is different from that of the averaged surface 
plane (z-axis in Figure [18]) . We introduce the polar coordinate with respect to the 

normal vector n of each facet as shown in Figure [181 Then the scattered radiance including 
the effect of rough surface is written as 



cos 9 



d cos 9 r , 



d cos 9' 



cos a; 



xS(9, 9', <f) - <p', 9 n , <j) n )R(u, u', f - Z% Q {B', <j)' 



(B3) 
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In the above, {9 n , <f) n } denotes the direction of the normal of the facet, and S is the effective 
fractional area of the wave facet associated with {6 n ,<j) n }. The integration over 9 n and 4> n 
results in setting oj — oj' due to the Dirac delta in Equation (1B2j) . Then Equation ( IBip 
reduces to 



outV u ) ' 



d cos 9' 
xS(9,9' 



<9(cos 6 n ,cj) n ) 



d(cos 



cos a; 



(B4) 



where 0* , and oj* are defined through the conditions of oj = oj' and £ — £' = 7r, or more 
explicitly, 



cos #* 
cos(2o;*) 



cos 9 + cos 6 1 ' 
2 cos 

cos # cos 9' + sin 6 1 sin 9' cos(</> — </>') 



as shown in Figure [THJ The Jacobian determinant in Equation (1B4|) is 

<9(cos 6> n ,0 n ) 



a(cos 



4 cos 



(B5) 
(B6) 



(B7) 



The effective fractional area S is expressed as 

E 



SS'd cos# r 



COS 61 



P{0 n , <j)n)T{8 , 6[; 4>, <f)'\0 n , <p n )dcos 6 n d<f) r 



(B8) 



where £ is the horizontal area, P is the density function of the wave slope, and T is the 
bidirectional shadowing fa ctor. The factor P(9 n , n ) is empirically given (see Figure 3 of 
Nakajima &: Tanakalll983l ) by 



1 



ira 2 cos 3 Q n 



exp 



tan 2 9 n 



a- 



P(0n) 



(B9) 



Here, a is the rms of slopes and a function of the wind velocity, uio, at 10 m height above 
the water surface taken to be 

a 2 = 0.00534w 10 . (BIO) 

Another factor, T(9,9';(/),(j)';9 n ,(j) n ), is obtained from an analytical fit to the results of 
a Monte Carlo simulation assuming one-dimensional random surface with Gaussian auto- 
correlation, which resulted in 



T(M';0,0'K, 



cos 9 
a sin 9' 
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1 + F(» +/>')' v 1 



F(v) = I 



exp(— t )dt 



(B14) 



where dz/dx is the slope of the facet, H{x) is the Heaviside step function, and a is given by 
Equation flBTUD . 



Combining all the expressions above, Equation (JB4J) is now written as: 

/out (e,^) = J dco&e' J d4>R{e,e\(j)-4>')i m {e'A'), (bis) 

mej, $-<!>') = ^—-p(e* n )G(v,v')r(u*,rh), (B16) 

cos 9 cos 9* n 

where 0* and u* are functions of 0, 0', and (0 — 0') through Equations (1B5j) and (1B6j) . 
Considering a parallel incident flux from the direction of {0o, 0o}, we have 

/ in (0,0) = FJ(cos6' - cos O )<5(0' - O ). (B17) 

Substituting Equation (IB 170 into Equation (1B15P and replacing {0, 0} with {6*i,0i}, we 
finally arrive at 

/ ou t(#i,0i) = F„i2(0 o , 0i, 0o-#i), (B18) 

^(00,01,01-00) = (B19) 

COS 01 COS 0* 

This is equivalent to Equation ( fl5|) in Section 12.2.41 : 

'—<*■ *> = W.coUcos^ ^ ">< W *- *>■ < B20) 



REFERENCES 

Arnold L., Gillet S., Lardire O., Riaud P., & Schneider J. 2002, A&A, 392, 231 

Baldridge, A. M., Hook, S. J., Grove, C. I., & Rivera, G. 2009, Remote Sens. Environ., 113, 
711 

Ballester, G. E., Sing, D. K. & Herbert, F. 2007, Nature, 445, 511 
Bennett, D. P., et al. 2008, ApJ, 684, 663 



-38 - 

Charbonneau, D., Brown, T. M., Noyes, R. W. & Gilliland, R. L. 2002, ApJ, 568, 377 
Cowan, N. B., & Agol, E. 2008, ApJ, 678, L129 

Cowan, N. B., Agol, E., Meadows, V. S. & Robinson, T. 2009, ApJ, 700, 915 

ESA SCI-A, 2007, Darwin Mission Summary Status Issue 2.0 (Noordwijk: ESA) 

Des Marais, D. J., et al. 2002, Astrobiology, 2, 153 

Ford, E., Seager, S. & Turner, E. L. 2001, Nature, 412, 885 

Frdlich, C, & Shaw, G. E. 1980, Appl. Opt., 19, 1773 

Goode, P. R., Qiu, J., Yurchyshyn, V., Hickey, J., Chu, M-C, Kolbe, E., Brown, C. T. & 
Koonin, S. E. 2001, J. Geophys. Res. Lett., 28, 1671 

Kasdin, N. J., et al., 2010, BAAS, 41, 287 

Easting, J. F., & Catling, D. 2003, ARA&A, 41, 429 

Easting, J. F., Whitmire, D. P. & Raynolds, R. T. 1993, Icarus, 101, 108 

Kiang, N. Y., Siefert, J., Govindjee, Blankenhship, R. E. 2007, Astrobiology, 7, 222 

Kiang, N. Y., et al., 2007, Astrobiology, 7, 252 

Lawson, P. R., et al. 2009, in Astro2010 — The Astronomy and Astrophys. Decadal Survey, 
53 

Levine, M., et al. 2009, "Terrestrial Planet Finder - Coronagraph (TPF-C) Flight Baseline 
Mission Concept". TarXiv: 091 1.32001 

Leger, A., Pirre, M., & Marceau, F. J. 1993, A&A, 277, 309 

Lucht, W., Schaaf, C. & Strahler, A. H. 2002, IEEE Trans. Geosci. Remote Sens., 38, 977 
Mayor, M., et al. 2009, A&A, 493, 639 

Menke, W. 1989, Geophysical Data Analysis: Discrete Inverse Theory (Revised ed.; New 
York: Academic) 

Montanes-Rodriguez, P., Palle, E., Goode, P. R. & Martin-Torres, F. J. 2006, ApJ, 651, 544 
Nakajima, T. & Tanaka, M. 1983, M. Quant. Spectrosc. Radiat. Transfer, 29, 521 



-39 - 

Oakley, P. H. H. & Cash, W. 2009, ApJ, 700, 1428 

Palle, E., Ford, E. B., Seager, S., Montanes-Rodriguez, P. & Vazquez, M. 2008, ApJ, 676, 
1319 

Press, W. H., Flannery, B. P., Teukolsky, S. A. & Vetterling, W. T. 1992, Numerical Recipes 
in C: The Art of Scientific Computing (Cambridge: Cambridge Univ. Press) 

Queloz, D., et al. 2009, A&A, 506, 303 

Salomonson, V. V., et al. 1989, IEEE Trans. Geosci. Remote Sensing, 27, 145 
Seager, S., Turner, E. L., Schafer, J. & Ford, E. B., 2005, Astrobiology, 5, 372 
Swain, M. R., Vasisht, G. & Tinetti, G. 2008, Nature, 452, 20 

Swain, M. R., Vasisht, G., Tinetti, G., Bouwman, J., Chen, P., Yung, Y., Deming, D., & 
Deroo, P. 2009, ApJ, 690, L114 

Tinetti, G., Meadows, V. S., Crisp, D., Fong, W., Fishbein, E., Turnbull, M. & Bibring, 
J.-P. 2006, Astrobiology, 6, 34 

Tinetti, G., Meadows, V. S., Crisp, D., Kiang, N. Y., Kahn, B. H., Fishbein, E., Velusamy, 
T. & Turnbull, M. 2006, Astrobiology, 6, 881 

Tinetti, G., et al. 2007, Nature, 448, 169 

Vidal-Madjar, A., Lecavelier des Etangs, A., Desert, J.-M., Ballester, G. E., Ferlet, R., 
Hcbrard, G. & Mayor, M. 2003, Nature, 422, 143 

Vidal-Madjar, A., et al. 2004, ApJ, 604, 69 

Wanner, W., Li, X. & Strahler, A. H. 1995, J. Geophys. Res. Lett., 100, 21077 

Williams, D. M. & Gaidos, E. 2008, Icarus, 195, 927 

Wolstencroft, R. D. & Raven, J. A. 2002, Icarus, 157, 535 

Woolf, N. J., Smith, P. S., Traub, W. A. & Jucks, K. W. 2002, A&A, 547, 430 

Young, A. T. 1980, Appl. Opt., 19, 3427 



This preprint was prepared with the A AS IATgX macros v5.2. 



